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Abstract 

In this paper we construct the Wilson short distance operator product expan- 
sion for the gluon, quark and ghost propagators in QCD, including operators 
of dimension two and three, namely, A 2 , m 2 , mA 2 , ipip and m 3 . We compute 
analytically the coefficient functions of these operators at three loops for all 
three propagators in the general covariant gauge. Our results, taken in the 
Landau gauge, should help to improve the accuracy of extracting the vac- 
uum expectation values of these operators from lattice simulation of the QCD 
propagators. 
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1 Introduction 



Two-point correlation functions of the fundamental fields of the QCD Lagrangian - that 
is gluon, ghost and quark propagators - are of direct importance in any perturbative 
treatment of QCD. Suffice it to say that the corresponding wave function renormalization 
constants are vital ingredients in calculations of the QCD /3-function and the quark mass 
anomalous dimensions (currently known at four-loop level [1-4]). Scheme-invariant ver- 
sions of these propagators are presently known in NNNLO (that is up to and including 
three loops) in arbitrary covariant gauge [5] including the Landau one, which is distin- 
guished from the point of view of lattice simulations. 

Purely perturbative treatment essentially assumes a weak-coupling regime. QCD prop- 
agators, especially the gluon and the quark ones, have been much under examination also 
beyond perturbation theory (that is in the strong-coupling regime). Here one should 
mention at least two broad directions, namely, the use of Schwinger-Dyson equations (for 
reviews see e.g., [6-8]) and non-perturbative computation on the lattice by Monte Carlo 
simulations. In what follows we will concentrate our discussion on the latter. 

It is expected — due to the asymptotic freedom — that the behavior of full QCD- 
propagators is to be governed at sufficiently large momentum transfers by perturbation 
theory and by the Operator Product Expansion (OPE) [9-11]. Thus, by comparing results 
of continuum perturbation theory calculations with those of lattice simulations one hopes 
to get a lot of information about the (renormalized) running coupling constant and quark 
masses as well as on condensates — Vacuum Expectation Values (VEV's) of composite 
operators — entering into OPE. The idea has been pursued in lattice simulations per- 
formed by various groups. (As for investigating condensates in lattice framework along 
these lines, see, e.g. Refs. [12-15] and also references therein for earlier results and for 
more lattice-specific information). 

While purely perturbative contributions to the QCD propagators have been computed 
in NNNLO, the corresponding (power suppressed) condensate contributions are usually 
known only at leading order or, at best, at next-to- leading order. To be specific, let us 
consider the gluon and ghost propagators in Landau gauge (for space-like momentum 
q 2 <0) 

sab 

D°(Q), A ab (q) = D h (Q) (1) 

with Q = a/— q 2 - The dressing functions D 9 and D h can be decomposed in terms of the 
appropriate OPE as follows 

D\Q 2 ) = DlWQ, a.) + £ CK *' Q d ; as) (Oi>, (2) 

i 

where ? stand for g or h, Dl(fi/Q,a s ) is the purely perturbative contribution, fi is the 

2 

renormalization scale and a s = ^ = is the quark-gluon coupling constant. The sum 
goes over all scalar operators with vacuum quantum numbers; di stands for the dimension 
of the operator Oi in mass units. 
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Assuming the case of massless QCD, the leading non-perturbative corrections in (2) 
should come from operators with the lowest possible mass dimension di = 2, namely (to be 
in agreement with the commonly used in lattice publication sign convention we effectively 
use below the euclidean scalar product in the definition of A 2 ) 

A 2 = -AlA a ^ and iC a C a , 

where A® is the gauge field , C{C) is the ghost (antighost) field. Within the class of covari- 
ant gauges 2 the coefficient function of the second operator is known to vanish identically 
in every OPE [16]. In what follows we will not consider this operator. 

The first operator, the gluon mass condensate^ 1 , does have nonzero coefficient functions 
C\ 2 already in the tree approximation, namely [16] (see, also [17-19]) 

C%=9 2 s ^ + 0(al), C h A2 =g 2 s ^ + 0(al). (3) 

From the phenomenological side, lattice simulations carried on by Boucaud et al in a 
series of publications [15,17,18,20-22] (see, also [23]) seem to demonstrate the existence 
of effects of order 1/Q 2 in both gluon and ghost dressing functions 4 . Moreover, numerical 
fits produce the results consistent with the OPE description of power suppressed 1/Q 2 
corrections to the gluon, ghost and quark propagators. This means, for instance, that one 
and the same value of the gluon mass condensate [15] 

g 2 (A 2 ) = 5.HJI GeV 2 (4) 

multiplied by the tree level CF's (3) together with purely perturbative contributions 
(known to three loops) describe the ghost and gluon dressing functions over the whole 
available momentum window 2 GeV < Q < 6 GeV. 

On the other hand, a study of a dressing function itself could, obviously, at best result 
in the determination of the product of the CF and the VEV of a composite operator 
(even if one assumes no contamination from operators of higher mass dimension). Thus, 
knowledge of higher order corrections to the coefficient functions of condensates is of some 
importance, at least for better understanding the results of lattice simulations. 

The quarks are massive. As a consequence the possible composite operators could 
contain powers of quark masses along with quantum fields. It is worthwhile to remember 
at this point that in "good" renormalization schemes like those based on the dimensional 
regularization [27-29] and minimal subtractions [30] the coefficient functions of any (short 



2 By a covariant gauge we mean the one generated by adding the term — 1 ^r L (d v Af l )(d u Af l ) to the gauge 
invariant Yang-Mills Lagrangian; the corresponding expression for the tree-level vector boson propagator 



reads 



; the choice of the Landau gauge corresponds to limit of £l — ► 0. 



g^v + (1 - ih)~ r 

3 We use this expression as the title of just gluon condensate is traditionally referred to the VEV of 
the operator G^G^ starting from the seminal works by the ITEP group [11]. 

4 The gluon mass condensate as well as the quark condensate also show up in the quark propagator 
[12,13,22,24-26], see Section 3 below. 



2 



distance) OPE obey the following important property 5 : their dependence on any parti- 
cle/field masses is polynomial. In particular, it means that any more complicated mass 
dependence of a correlator will be "hidden" in the corresponding VEV of composite oper- 
ators. It also means that if one allows, as we do, mass factors to be used in constructing 
composite operators, then their coefficient functions become totally mass-independent by 
definition. Limiting ourselves to operators with mass dimensions not higher than three 
we arrive to the following list of operators which could appear in OPE for the QCD 
propagators: 

A 2 = A^A^, m 2 , m 3 , mA 2 , (5) 

where m is a quark mass, and we have assumed QCD with rif — ni + 1 total number of 
quark flavors, one of those having a mass m, while all others are strictly massless 6 . 

The aim of the research we are going to present was to compute the higher order 
contributions (up to and including three loops ) to coefficient functions of operators (5) 
appearing in the OPE of QCD propagators. The structure of the paper is as follows. In 
the next two sections we describe our results for OPE of QCD propagators. In the fourth 
section we consider the RG evolution equations for propagators and operators under 
consideration and construct the scale and scheme invariant combinations of operators 
and coefficient functions. Due to their scheme independence the latter should be most 
convenient for comparisons with the results of lattice calculations. Then we briefly discuss 
(in Section 5) some technical details of the calculations as well as software/hardware tools 
employed. Finally, a short summary of our findings is given in the Conclusion (Section 6). 

We finish the introduction by adding that in recent years, starting from works [33,34], 
the condensates of mass dimension two, especially the gluon mass condensate, have been 
intensively studied in view of better understanding of confinement in Yang-Mills theories 
and QCD. (For example, see recent works [35-43] and references therein). Unfortunately, 
any discussion of these developments is beyond the scope of the present paper. 

2 OPE for the gluon and ghost propagators 

On dimensional grounds, the operators of dimensions three do not contribute to the OPE 
for the gluon and ghost propagators. The remaining coefficient functions C^ 2 and C\ 2 

5 To our knowledge it was first established in [31]; see also [32]. 

6 Later, in Appendix A, we generalize our results for the case of arbitrary many massive quarks. 
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Here and everywhere in the paper the renormalization is carried out in the MS-scheme, 
rif is the total number of quark flavours, I^q = In £5, m = m(«) and a s = are 
the running quark mass and quark-gluon coupling constant respectively. In addition, the 
irrational constants £3 = 1.2020569, (4 = 1.0823232, £5 = 1.0369277 appear. In numerical 
form Eqs. (6)- (9) read 

C g m2 =a s [I + a s (17.7614 - 0.555556 n f ) 

+ a*(311.276 - 25.8972 n f + 0.231481 n))] , (10) 

C% =- n 2 a s [I + a s (8. 17708 - 0.611111 n f ) 

+ 0^(86.9600 - 13.6255n/ + 0.373457 n)) 

+ a^(1068.69 - 240.803n/ + 14.1703nJ - 0.1911874)] , (11) 

C* a = - -a 2 s [1 + a s (33.2934 - 0.958333 n/)l , (12) 
8 

C% =ln 2 a s [l + 3. 75000a, 
8 

+ ^(30.0240- 1.4609471/) 

+ a*(339.141 - 38.2844 n f + 0.655878 n))] , (13) 
where we have set ji = Q. 



3 OPE for the quark propagator 

The quark propagator of a quark field ip q with mass m is expressed in terms of the 
corresponding dressing functions as follows: 

i J dse*<TMs$(0)]> = ±V(Q) + ^ (14) 
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The OPE expansions for the dressing functions (up to operators of dimension three) are 
V{Q) = Vo(ti/Q,a s ) + m q 2 m + A (A), (15) 



S{Q) = s „ ( ,/ e ,a, )m + m 3 + + 

(16) 

where m is the quark mass of the quark associated with the quark field ip. 

The purely perturbative contributions Vq and So have been already discussed at three- 
loop level in [44], the files with results in computer readable form can be downloaded 
from http : / / www-ttp . particle . uni-karlsruhe . de/Progdata/ ttp99/ ttp99-43. The 
remaining coefficient functions are listed below: 
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Their numerical form (with (j, — Q) reads: 

C q m2 = - [1 + 2.66667a s + ^(21.7015 - 0.840278 n/) 

+ a^(221.404 - 19.7294 n f + 0.173079 n))] (22) 

7T 2 

C% =- —a a [l + 0.750000a, + ^(2.39159 - 0.178385 71/) 

+ a;j(20.2621 - 2.20883 n f + 0.0377513 n))] (23) 
C q mi = - [1 + 4.00000a s + ^(36.1257 - 1.25000 n/) 

+ aij(397.959 - 34.3465 n f + 0.397359 n})] (24) 

257T 2 

C Ia- =-4g- a 2 + a.(H-0225 - 0.414444 n/) 

+ ^(135.998 - 12.805972/ + 0.185864 n))] (25) 

471-2 

= - — a s [l + a s (6.18750 - 0.277778 72/) 
+ a 2 (52.9495 - 6.52486 n f + 0.0771605 n)) 

+ a;j(564.828 - 112.308 n f + 4.70773 n) - 0.0214335 n))] (26) 

4 Renormalization Group Improvement 
4.1 Anomalous dimensions 

We limit ourselves to the three-loop level. Let us start from the operators of dimension 
two. The corresponding matrix of anomalous dimensions is defined by the following matrix 
equation 

d ( A* \ ( 7A2 7AW \ ( A 2 



where the differentiation on the lhs is carried out with fixed bare coupling and quark 
masses. The quark mass anomalous dimension is known since long [45,46] and the anoma- 
lous dimension of A 2 in Landau gauge was found in [47] to be 
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The non-diagonal three-loop anomalous dimension 7^2 , m 2 reads 

'971 



7a 2 ,m 2 



16tt 2 



24 + a. 



-471/ + 36C3 



(29) 



Life is easier with operators of dimension three. First, the anomalous dimensions of the 
pair m A 2 and m 3 are, obviously, additively related to those considered above, namely: 



7m A 2 = 7m + 7,42 , 7m3 = 37„ 



(30) 



Second, in the process of renormalization the quark condensate could mix only with the 
unit operator (times a quark mass cubed): 



d 

dy 



m 3 



37 m 



m 3 



(31) 



The fact that 

7^ = ~7m (32) 

is well-known from text-books. The non-diagonal part of the mixing was investigated in 
detail a long time ago [48,49]. It is naturally expressed in terms of the so-called vacuum 
anomalous dimension, 7q as follows [50]: 



with 



7o d 



A* J^2^ = ~7m ^ - 4m %(a s ), 



, 4 /313 5 2 A . ., 

i + ^+l^-iYv-gCsK 



16tt 2 



(33) 



4.2 Scheme-independent correlators and operators 

In general a Green function G depends on the renormalization prescription (scheme) and 
the choice of the artificial scale fi. It is, however, well-known how to define a variant G of G 
which is invariant under changes of the renormalization scheme and /i. The corresponding 
renormalization group equation (RGE) 



fi 2 —G = 



d/i 2 



has the formal solution 



G — G(a s ,fj)/f(a s ), f(a s ) = exp 



f dc 
J a 



da s jc 



(34) 



(35) 



The OPE of a (suitable) scale-invariant Green function can be rewritten in terms of scale- 
invariant operators Oi and Wilson coefficients C i: which again obey RGEs of the form 
(34). 
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Here, we consider the OPEs of the scheme-independent dressing functions in the limit 
of massless quarks. In this limit — aside from the perturbative contributions — only 
the operators A 2 and ipip contribute. For the operator A 2 and its coefficient functions we 



obtain (in the three-loop approximation) 

_ 35 

A 2 \ nj=Q = a s 44 (1 + 0.0693440 a s + 0.0240863 a 2 s + 0.405494 a 3 s ) A 2 , (36) 

A 2 \ = a~ s ^ (1 + 0.0654912 a s + 0.0933818 a 2 s + 0.508904 a 3 ) A 2 , (37) 

A 2 ! =aj*(l + 0.0538194a s + 0. 136131a;? + 0.570436 a 3 s ) A 2 , (38) 

C%\ n=0 =\ -k 2 cl}~ a (a s + 8.24643 a 2 s + 87.5512 a 3 s + 1075.32 a 4 ) , (39) 

C%\ nf=2 = - vr 2 a s ^ (a s + 7.02035 a 2 s + 61.7518 a 3 s + 647.400 a 4 ) , (40) 

^ |„ = 3 = r 71-2 a * ( a * + 6.39757 a 2 + 49.9224 a 3 s + 472.744 a 4 ) , (41) 

C% \ nf=Q = - v 2 af (a s + 3.61131 a 2 s + 29.4702 a 3 s + 334.048 a*) , (42) 

C* 2 \ n =2 = - 7r 2 aF (a s + 3.61902 a 2 + 26.4370 a^ + 260.012 a 4 ) , (43) 

C% 1^ =3 = - 7T 2 v^I(a s + 3.64236 a 2 + 24.9740 + 225.345 a 4 ) , (44) 

■jt-2 35 

C %\ n =o = ~ -^ a * 4 ( a s + 0.173080a 2 + 1.19104 a 3 s + 15.9766 a 4 ) , (45) 

C%\ nf=2 = - ya s ^ (a, + 0.175888 a 2 + 0.859625 a 3 s + 12.6022 a*) , (46) 

7T 2 3 

C* 2 | n =s = a /(a s + 0.186921 a 2 + 0.689610 a^ + 11.0348 a 4 s ) (47) 

by inserting the corresponding anomalous dimensions into Eq. (35). The scale- and 
scheme-independent versions of ifjifj and its coefficient function in the OPE of the quark 
propagator read 

U| = (1 + 0.687328 a s + 1.51211 a 2 + 4.05787af) ^ , (48) 

_ =Os^ (1 + 0.805985 a s + 1.40095 a 2 + 2.72916 a 3 .)^, (49) 

Hj\ nf=3 =a~H l + 0.895062 a, + 1.37143 a 2 + 1.95168 a^)^, (50) 

CJJ B/=0 = - —a/ 1 (a. + 4.99260 a 2 + 44.0815 a* + 489.206 aj) , (51) 

47T 2 1? 

C;-J n/=2 = - — a » (a s + 4.31734a 2 + 31.7744a* + 298.894 a*) , (52) 

47T 2 i 

CJ | n =3 = a/ (a, + 3.94985 a 2 + 25.7972 a 3 s + 217.583 a*) . (53) 

f 3 
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Another useful scheme-invariant object is the so-called "effective quark mass" m P {Q) 
which is defined as follows [10]. 



i / dxe iqx (T[^(x)ij(0)}) = 



1 B{q) S{Q) , KA . 

b^it€ mm = m = mr (54) 

where we have used Eq. (54) to express mp(Q) in terms of the dressing functions V(Q) 
and S(Q). 

In the chiral limit the leading contribution to mp(Q) comes from the quark condensate; 
in explicit form we get 

rn P (Q) = C^V -^. (55) 



Using the results of Section 3 we arrive at: 
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m P (Q)\ n 2 = - 
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(56) 
(57) 
(58) 
(59) 
(60) 
(61) 
(62) 
(63) 



where we have set the renormalization scale ji = Q. Numerically these equations read: 



mp(Q)| n/=0 
m P (Q)\ nf=1 
mp(Q)\ n}=2 
m P (Q)\ nj=3 



47T 2 

3 

47T 2 

3 

47T 2 

3 

47T 2 



a s (l + 6.1875 a s + 54.541 a 2 + 597.957 a 3 ) 



Q 2 



a s (l + 5.90972 a s + 47.9475 a 2 + 485.281 a 3 s ) ^ ' 



Q 2 ' 
Q 2 ' 



a,(l + 5.63194 a s + 41.5083 a 2 + 382.176 ) 
a s (l + 5.35417 a s + 35.2234 a 2 + 288.511 a 3 s ) 



s Q 2 



(64) 
(65) 
(66) 
(67) 
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In analogy to the effective quark mass one can also define effective masses for gluon 
and ghost fields, which are induced by the gluon mass condensate (A 2 ). An explicit 
formula can be derived by considering the ghost propagator (or the gluon propagator) in 
the chiral limit: 



A ab (q) 



6 ab 

Q 2 



D h (Q) + c h AQ) 



Q 2 



5 



ab 



Q 2 -C h A2 (Q)/D^Q)(A 2 ) 



D o(Q) • (68) 



The effective masses are then given by 

m 2 (Q) 



&AQ) 



{A 2 



D l(Q) 

where ? stands for g or h. The analytic results for /i = Q and rif — 1, 2, 3 are 



(69) 



rn 



g\n f =0 



m 



g\ n f= l 



m 



TT 2 a s 



8 

+ a 3 s 

+ aj^- 
35721 

+ 



l + a & 



197 



32 



2 /480587 243 
+ C H^2T6" + T28 C3 



/ 520248245 



V 884736 4096 



82215 , 243 , 331077 , 



4096 

38766561211 25948995 „ 



8192 



7077888 



131072 



57028439 729 2 
32768 ^ 4 + 131072 ^ 3 + 2048^ 3 
3 9 T /559 

( 11737602763 74855 



(A 2 



7r a s 

8 



Z3836375 195 
1 82944 + 128^ 3 



V 23887872 



4096 



285 , 7127207 , 



4096 



221184 



2182125807293 1465676035 
+ 



573308928 



38475 , 12243145349 , 
C 4 + G 



10616832 
8075 



<5 



32768 
3 



31850496 
527 
~96~ 

/9610131071 67495 



-ir 2 a 



1 + a 6 



18432 
3368459 



(A 2 ) 
147 

82944 + 12lr 3 



V 23887872 4096 



(5 



813 „ 5402375 „ 



4096 



477599196341 916103635 A 



191102976 



99999 , 3506053415 , 
^4 + An „ nonn C 3 



10616832 
24119 



32768 



10616832 



18432 



221184 



(A 2 



(70) 



(71) 



(72) 
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3 2 r /165\ ,/ 108205 99 . \ 

7i 2 a s 1 + a J + ail 1 C3 

8 3 1 V 32 / V 3072 128 S 7 

g/31529153 _ 60135 1341 139429 
+ ° s V 98304 4096 ~ 4091) + ~8192~ 

J 10684895233 16783385 
+ a \ 7077888 + 393216 
148851 322973567 4619 2 \1 2 
32768 ^ 4 + 1179648 ^ 3 ~ 2048*7 J ^ ' ' 



2 
-7r a s 
8 



2 / 20997 
1 + 3a ' + a HT024 



2 / zuyy 1 351 \ 
5 V 1024 + 256^7 
3/ 4345483 53235 243 24945 
+ a ° V 18432 ~ 4096 ^ 5 + 4096^ 4 + 1024 ^ 3 

+ 



, 4uyo 4uyo iuz 

4 / 1004316139 222885 5103 4937751 3645 2 Y 
° s I 1572864 + 16384 ^ 5 ~ 16384 ^ 4 + 131072 ^ 3 + 8192 *y_ 



3 2 [ 2 / 19881 351 \ 

-&* a 4 1 + 3as + as lW + 256 C3 J 

g/34790023 53235 285 69767 \ 
+ ° s V 165888 ~ 4096 ^ 5 ~ 4096^ 4 + 3072 ^) 

J 8001655387 222885 5985 4761971 3645 
+ ° s V 14155776 + 16384 ^ 5 + 16384 s4 + 131072 ^ 3 + 8192' 



7r a s 

8 



1 + 3 a s + a„ 



18765 351 



1024 



+ 



73( - ' • 53235 813 



V 165888 4096 



<5- 



4096 



256 

Ci 







64763 



3072 



Cs 



+ a* s \ - 



7006421587 222885 
+ 



14155776 



16384 



C5 + 



17073 
16384 



C4 + 



4571855 
131072 



, 3645 , 

C3 H C 

S3 8192 s ' 



J 2 
7T a 

8 



351 N 
256 s3 



' 2 / 17649 
1 + 3 a s + aj 

3/ 329079 53235 1341 19941 

+ ° s V 2 048 ' 4096 ^ 5 ~ 4096^ 4 + 1024 ' 

J 224190513 222885 , 28161 
+ a 4 J + Cs + tz^ttCa - 



' 4096" 1024 "7 
222885 28161 4367403 3645 s 
524288 + 16384 ^ 5 + 16384 ^ 4 + 131072 ^ 3 + 8192 s3 



( 
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In numerical form we obtain 



2| 

m g \n r - 


=0 — 


3 *n 

7r a s 

8 


l + o. 15625 a s + 54.4291 a s + 615.858 a s — 4747. 15 a s \ {A ) , 


(78) 


2| 

m g\n r - 


=1 = 


3 2 

7r a s 

8 


1 i cr ooono i ao noon 2 i ri i n^i 3 oono on 4] / /t2\ 

1 + 5.82292 a s + 48.0839 a s + 511.071 a s — 3202.89 a s \ {A ) , 


(79) 


2| 


=2 = 


3 2 

7r a s 

8 


1 + 5.48958 a s + 41.9917 a s + 414.36 a s — 2017.93 a s \ {A ) , 


(80) 


2 1 


=3 — 


3 2 
8 


1 + 5.15625 a s + 36.1527 a s + 325.612a^ - 1144.42 a s \{A ) , 


(81) 




=0 = 


3 2 

7r a s 

8 


1 + 3a s + 22.153 a\ + 251.628 a\ - 578.831 a*] (A 2 ) , 


(82) 


m h\n f -- 


=1 — 


3 2 

7r a s 

8 


1 + 3 a s + 21.0632 + 223.467 a 3 s - 506.441 a*] (A 2 ) , 


(83) 


m h\n f -- 


=2 = 


3 2 

7r a s 

8 


1 + 3 a s + 19.9733 + 196.345 a 3 s - 437.146 a*] (A 2 ) , 


(84) 


m h\n r - 


=3 — 


3 2 
7r a s 

8 


1 + 3 a s + 18.8835 + 170.26 a 3 s - 370.947 a*] (A 2 ) . 


(85) 



5 Theoretical and hardware tools 

In our work we have heavily used the software packages QGRAF [51], EXP [52] and a 
modern version of MINCER [53,54] written in the algebraic computer language FORM 3 
[55] for the generation and calculation of the required diagrams. 

The calculations have have been performed in a general covariant gauge 
and for the gauge group SU(n). The total number of diagrams contributing 
to different channels (according to QGRAF) are displayed in Table 1. The 
files with all results (in a computer readable form) can be downloaded from 
http : //www-ttp . particle . uni-karlsruhe . de/Progdata/ttp09/ ttp09-40/ 

Any coefficient function of any operator entering into the OPE of two local operators 
can be expressed in terms of massless propagators. The reduction to massless propagators 
is conveniently done with the well-known method of projectors [56,57]. 

Finally, the package MINCER is able to compute very effectively massless propagators 
up to (and including) three loop level. 

6 Conclusion 

We have computed the coefficient functions of the operators of dimension two and three 
in the OPE for the gluon, ghost and quark propagators. The higher order corrections 
are essential as one could see by inspecting eqs. (36) -(53). They are most important 
in two cases: gluon condensate contributions to ghost and gluon propagators as well as 
the quark condensate one to the quark propagator. In general the terms proportional to 
rif tend to significantly stabilize the perturbative series by decreasing the value of higher 
order terms (cmp. e.g. eq.(39) and (42)). 
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tree one loop two loops three loops 




1 5 59 1148 
6 222 7407 264399 
11 9 148 

2 23 595 19419 
11 9 148 
2 23 657 23251 
1 11 234 6641 



Table 1: Number of diagrams contributing to the various coefficient functions 

Specific numerical analysis should be made with specific lattice data. Still, we observe 
that the higher order corrections to the coeffcient functions display relatively good (ap- 
parent) convergency pattern with basically positive coefficients which, presumably, should 
lead to a noticeable decrease of the value of the A 2 condensate once the lattice data are 
reanalyzed with an account of newly computed terms in the corresponding OPE. 

Note that such dependence of the numerical value of the gluon mass condensate on the 
number of perturbative terms accounted in the corresponding OPE gives an extra support 
to the the hypothesis of duality between perturbative and non-perturbative contributions 
(see a very recent work [58] and references therein). 

Finally, we hope that our results will be of use for better understanding of the present 
and future data coming from lattice simulations of QCD propagators. 
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A Several massive quarks 

Until now, all OPEs in this work have been formulated for the case of (at most) one 
massive quark with mass m. From a physical point of view, this is a valid approximation: 
corrections from u and d masses are in general negligible, while all other quark masses 
show a strong hierarchy. 

Nevertheless, it is also possible to generalize our results to the case of n/ massive 
quarks with masses m l5 m 2 , ... ,m nf . The generalization of the OPEs of the gluon and 
ghost propagators is straightforward: We replace m 2 by the sum over m 2 : 



D\Q 2 ) = D 7 MQ,a s ) + 



C ? A 2{y/Q,a s ) 

Q 2 



(a 2 ) + E 



n 



^C 7 m2 (fi/Q,a s ) 




(86) 



i=i 



Q 2 
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where C^ 2 is the same as C 7 m2 from Eqs. (6) and (8). 

The case of the quark propagator is a bit more complicated. Without loss of generality 
we assume that the external quark has the mass m\. The correspondingly generalized 
OPE for the dressing functions V(Q) and S(Q) defined in Eq. (14) then read 



V{Q) = Voifi/Q, a,) + 2^ -7y2 + Tyz 



Q 2 1 Q 2 



(87) 



and 



^C q m2mi (fi/Q,a s ) 

S(Q) =S (fi/Q,a s )m 1 + ^ — l - — — m.mi 

i=i V 

+ ) + q 2 (W ■ 

Note that for % ^ 1 only even powers of rrij may appear in the OPEs. 

The Wilson Coefficients C q 2 , C" 3 3 and C 9 , 2 are identical to C q 2 , 3 and C" 3 , 2 as 
given in Eqs. (17), (19) and (20). For % ^ 1, the anomalous dimension matrix of (31) has 
to be extended to account for the mixing of ipifj with m\m\. The corresponding anomalous 
dimension 7^, m ? mi can be found in Ref . [50] . The results for the new coefficient functions 
are (note that below i ^ 1!) 

C^ 2 =a 2 (--\l q 



"s 



6 2 

3 / 63659 2243 129 2 5 3 



3456 192 64 v 2 



1. 3 A r 37 13 T 1 r2 . 



3\ , ^17143 319 



C* =a 2 A-\+all^-^-^L 



m^-us y2)^ a ° v ^76~~^2 Q 

-5Cs + ^Ca - \n f + \n f L Q ] . (90) 
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